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Abstract 

Systems consisting of cold interacting bosons show interesting collective phenomena such as Bosc-Einstcin condensation or 
superfluidity and are currently studied in condensed matter and atomic physics. Of particular interest are nonidcal bosons 
which exhibit strong correlation and spatial localization effects. Here we analyse the ground-state of a two-dimensional Bosc 
system with a Hartree-Fock type approximation that was first introduced by Romanovsky et al. [Phys. Rev. Lett. 93, 230405 
(2004)]. We apply this method to a one dimensional system of charged bosons and analyze the behaviour at strong coupling. 
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1. Introduction 

With the first experimental realisation of a Bose- 
Einstein condensate in 1995 [1] , the theory of ultracold 
bosonic systems became a field of exceptional interest. 
In the first experiments, the atoms were weakly inter- 
acting, and for such systems the Gross-Pitaevskii mean 
field approximation (GP) is an adequate approach for 
the analysis of the condensate [2]. However, recent ex- 
periments succeeded in strongly increasing the interac- 
tion strengths of the investigated atomic systems [3]. 
With this tuning of the coupling, it was possible, to 
produce a phase transition from the superfluid to the 
Mott insulator phase [4]. Also, for charged bosons in 
traps crystallization and inhomogeneous distribution 
of the superfluid density has been predicted [5]. For a 
theoretical description of these effects an approxima- 
tion beyond GP is required [6,7,8]. In this paper we 
analyse an approximation method for bosonic systems 
which is based on a Hartree-Fock type factorization 
and is expected to be valid for a wide range of interac- 



tion strength. 

In the many-body theory, Hartree-Fock is a standard 
method to analyse systems of few to many fermions. 
The related equations can be derived in several differ- 
ent ways. One is obtained by expanding the reduced 
two-particle density matrix in terms of one-particle 
density matrices in order to obtain an effective single- 
particle hamiltonian. The resulting solutions to this ap- 
proach are many-particle states consisting of one sin- 
gle anti-symmetrised product-state. Another common 
way proceeds in the opposite direction by the assump- 
tion, that the many-particle state is given by a single 
anti-symmetrised product state (Slater determinant). 
By applying the Ritzian principle to this ansatz, one 
obtains the Hartree-Fock equations. 

In contrast to Fermi systems, for bosonic particles 
the two approaches lead to different systems of equa- 
tions. Even the ansatz assuming a single Slater perma- 
nent can be implemented in various ways differing in 
the choice of supplementary conditions [6,7]. 

The best known approximation for an ultracold non- 
ideal dilute Bose gas is the Gross-Pitaevskii approxi- 
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mation in which it is assumed, that all particles remain 
in the same one-particle orbital. Thus the GP equa- 
tion describes only the condensate fraction of the sys- 
tem without any interaction with the remainig gas. By 
definition the GP approximation cannot describe phe- 
nomena such as condensate depletion and fragmenta- 
tion or Mott-insulator phase transitions of cold atomic 
gases on a lattice. The results become even worse if the 
interactions increase or if systems with small particle 
numbers are considered. As we will see later, this ap- 
proximation cannot describe a number of fundamental 
properties, such as the localisation of the particles at 
high coupling strengths. 

Another common numerical method to deal with in- 
teracting systems is the Configuration Interaction (CI) 
in wich no approximation, except the limitation of the 
chosen basis set, is made, e.g. [9]. With nowadays com- 
puters and highly optimised programs CI can be ap- 
plied to systems with up to 10 particles. 

In this paper we want to analyse an approximation 
with substantially less numerical effort than CI, that 
provides high quality numerical results for strongly cor- 
related bosons in traps. This approximation, the Un- 
restricted bosonic Hartree-Fock method (UBHF), was 
first introduced by Romanovsky et al. [10,7,11] and 
is less restricted than the related ansatz proposed by 
Cederbaum [6] . While Romanovsky et al. used an ex- 
plicit analytical ansatz for the single-particle orbitals 
(discplaced Gaussians) we develop a completely gen- 
eral scheme without any such restriction. This has the 
advantage that our method is applicable to any inter- 
acting Bose system. 

Furthermore, our goal is to compare the UBHF re- 
sults with CI results to obtain a quantitative conclu- 
sion about the accuracy of this method. To this end 
we perform UBHF and CI calculations for N = 2 ... 6 
particles in a harmonic trap using the same basis sets 
for both cases. Our comparisons concentrates on sys- 
tems with a coupling parameter in the range from zero 
to five and shows the excellent quality of this approxi- 
mation. 



2. Unrestricted bosonic Hartree-Fock ansatz 
(UBHF) 

The UBHF method is derived from the following 
ansatz for the many-particle state 

|$) = \12...N) := — L= (|)kW>. (1) 

which means that every particle remains in a certain or- 
bital and the underlying many-particle state is a sym- 
metrised product-state. Depending on the imposed ad- 
ditional restrictions, this ansatz contains, as limiting 
cases, well-known approximations. In particular, the 
GP approximation is obtained by the additional re- 
striction, that all orbitals are identical [8], thus the 
many-particle state is assumed to be totally Bose con- 
densed. A more general approximation would be ob- 
tained by the assumption, that two orbitals are either 
equal or orthogonal to each other. This kind of ansatz 
has been introduced and analysed extensivly by Ceder- 
baum et al. [6,12,13,14]. 

In contrast, the UBHF approximation is the most 
general case of this Ansatz where no further restriction 
to the underlying one-particle orbitals |1) , . . . , \N) is 
imposed. We only require the many-particle state to 
be normalised, so the total energy is given by 

e = (<z>\h\<s>) . (2) 

Thus minimising E leads to the functional 

E(\l),..., \N) , S) = <*|ff|S) - <?«$]$> - 1), (3) 

where § is a Lagrange multiplyer for the normalisation 
of |$). The needed equations to determine the one- 
particle orbitals are obtained by applying the Ritzian 
principle to this Ansatz via the searched orbitals. This 
ansatz has been proposed and implemented for the first 
time by Romanovsky et al. [11,10], however, they used 
the additional assumption that every orbital is a dis- 
placed Gaussian. In the present paper, this restriction 
will be dropped. 

The expectation value of the energy of a normalised 
single permanent many particle state with interacting 
particles is given by [15] 

1 — N 

way*) = -±= Edl^wx'i^')) 

+ ( S k( S )}«'fei, 7 r(fc)vr(!))- (4) 
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The norm of such a productstate is proportional to 
the permanent of the Gramian matrix of the orbitals 
|1>,...,|JV>: 

1 N 



It is advantageous to stay in the abstract representa- 
tion of the orbitals to perform the variational deriva- 
tive. By applying the rules given in Appendix A one 
obtains the following equation, by differentiating (3) 

]T I P n h |ir(n)> + p «i <1|&M0> K»)> 



= S E P ™ l""^)) ' fol ' ^ n ' 

where we introduced the subscripted entity with a vari- 
ational number of subscripts 

iv.<»=-p(T)ii...<»= n < s M s )>> w 

and the doubly subscripted Hartree operator 



Jij = J dxdy tpt (x)w(x, y)tpj (x) \y) {y\ . (8) 

By multiplying equation (6) from the right with (n[, 
one obtains a closed expression for the Lagrange mul- 
tiplyer 

mm) 



s 



($1$) 



0) 



i.e $ is the total energy of the System. The key equation 
(6) can be solved with a multidimensional minimisation 
routine without making any assumption on the explicit 
analytical form of the orbitals. 
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Fig. 1. Comparison of the energies obtained with UBHF, 
GP and CI for different coupling parameters A and for 2,3 
and 4 particles. 



x is given in oscillator units xo = ^/h/muj, and A is the 
ratio of the Coulomb and confinement energy 



A : 



(11) 



Thus A = corresponds to an ideal system. The shield- 
ing parameter re is needed in ID to make the two- 
particle integrals convergent. For the limit re — > 
the system becomes fermionised [16]. For all calcula- 
tions in this paper, we choose re = 0.1. As basis for the 
single-particle orbitals we chose the eigenfunctions of 
the ideal system with a total number of basis states of 
rib = 15. We will compare the results with the exact 
solution obtained with CI [9] . The CI calculations are 
done in exactly the same basis, thus the differences in 
the results arise exclusively from the UBHF ansatz. 



3.1. Total energy 



3. Numerical results 



We applied the UBHF approximation to few inter- 
acting bosons in a one-dimensional harmonic trap de- 
scribed by the hamiltonian 




which has been made dimensionless by using standard 
oscillator length and energy scales: the spatial variable 



In Fig. 1 the energies obtained by both methods are 
compared. Interestingly, the difference between the ex- 
act method and the UBHF approximation becomes a 
constant shift for high interaction values, A > 2. This 
shift still depends on the considered particle numbers 
N. As one can see in 1, the shift grows with N. Note 
that for GP the energy diverges rapidly from the exact 
result already for A > 0.5. Thus Fig. 1 is a convincing 
evidence for the good quality of the UBHF approxima- 
tion. 
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Fig. 2. Particle density of a system with N = 2 and A = 0.9 
for the three different approximation methods. The two 
orbitals obtained with UBHF are also given. 




Fig. 3. Same as Fig.2, but for A=1.5 
3.2. Localisation of the orbitals for X 2> 1 

The UBHF approximation scheme offers the possi- 
bility to analyse the delocalisation of the interacting 
particles in a special way. As can be seen in Figs. 2, 
3 and 4, the overlap of the orbitals vanishes with in- 
creasing A. It is crucial to notice, that the density from 
UBHF shows a localisation of the particles which, by 
definition, is missing in the GP model. Finally, the 
other two curves in Figs. 2 and 3 which resemble Gaus- 
sians are the two orbitals obtained by UBHF. It is in- 
teresting to see that the localisation of the particles 
emerges already in each of the orbitals. This can be 
understood as a precursor of the classical strongly cor- 
related limit of the system where the particles form a 
fully localized crystal-like arrangement [17,18]. 

This trend is typical for the present system of 



Coulomb interacting trapped bosons. As another ex- 
ample, in Fig. 4 we show the density of six bosons 
for two values of the coupling parameter. Again, with 
increasing A the overlap of the orbitals decreases and 
the system, as a whole expands. 

3.3. Delocalisation and non-orthogonality of the 
orbitals for small A 

There are several possibilities to analyse the overlap 
of orbitals in a quantitative way. One that works for 
all particle numbers is to consider the Gramian deter- 
minant of the given set of orbitals, 



G(\l): 



,\N))= J2 sigaM n <>!*(*)> 

<1|1> <1|2> ... (1\N) 
(2|1) <2|2) ... (2\N) 

(N\l) (N\2) ... (N\N) 



(12) 



This entity is always positive and approaches the value 
1, if the orbitals form an orthonormal set (this is most 
easily seen for the case of two orbitals). Geometrically, 
the Gramian determinant is the square of the volume of 
the parallelepiped spanned by the vectors |1) , . . . , \N) 
[21]. In Fig. 5 the A-dependency of G is shown. Obvi- 
ously, in the limit of large A the orbitals become pair- 
wise orthogonal. With increasing particle number this 
limit is reached for larger values of the coupling param- 
eter. 

Let us now consider the opposite limit of small cou- 
pling. As Fig. 5 shows, for all N, the determinant G 
monotonically decreases when A is reduced. Since G is 
a determinant it is obvious that it vanishes, if one or 
more orbitals are collinear. In the present case, also the 
opposite is true: vanishing of G is an indication of the 
orbitals becoming collinear. This is observed for A — > 
0. Due to the symmetry of the system not only two or- 
bitals but all orbitals are becoming collinear, i.e. they 
are identical. But this is just the limit of an ideal Bose 
gas where we expect that all particles Bose condense 
in the ground state. It is a remarkable property of the 
present UBHF ansatz which does not impose any re- 
strictions on the orbitals that it yields the correct BEC 
limit of identical Hartree Fock orbitals for all particles. 
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Fig. 4. Connection of the obtained orbitals (upper two fig- 
ures) and the resulting total particle density of a system 
with N = 6 particles and A = 0.5 and 1.0. 
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Fig. 5. Gramian determinant of the set of orbitals for dif- 
ferent particle numbers versus A 

4. Discussion 

In this paper we have studied a Hartree-Fock ap- 
proximation for interacting bosons in a trap which has 
proposed by Romanovsky et al. In contrast to their 
work we have not used any assumption on the explicit 
form of the orbitals but obtained them selfconsistently. 
We have tested the UBHF scheme in detail by compar- 
ing the resulting energies and densities with the ones 
obtained with CI. Due to the identical choice of the 
basis set for both approximations, the differences in 
the results arise exclusively due to the approximation 
made in UBHF. The comparison shows that this is a 
very accurate model. In contrast, it is confirmed that 
the GP approximation is not applicable to nonideal 



charged bosons in traps when the dimensionless inter- 
action strengths A exceeds 0.5. For its resulting ener- 
gies are substanially higher than the exact ones and it 
cannot reproduce important qualitative features of the 
system such as the localisation of the total particle den- 
sity. In the limit of vanishing A the UBHF approxima- 
tion yields single particle orbitals which are collinear - 
in other words, we automatically recover the GP model. 

We mention that the computer time required to solve 
the equations (6), within our current implementation, 
has an unfavorable dependence on the particle number 
N scaling as N ■ N\. This arises from performing the 
sum over all permutations in (6). However, it is possi- 
ble to reduce the calculation of the entities that turn 
up in this equation to the calculation of permanents. 
By using the Ryser algorythm, the JV-dependency of 
the complexity becomes of order TV ■ 2 N [19]. This is 
still a rapidly growing calculation time, but it is much 
better than the first version. The calculation time also 
depends on the size of the chosen basis set «,(,. For both 
implementations it is of order nf. This arises from the 
calculation of the two-particle-integrals in each step 
[20]. 

Presently we are working on further optimizing the 
implementation of this ansatz in order to extend it to 
systems with higher particle numbers and higher di- 
mensions. Furthermore, the development of an exten- 
sion to the time dependent regime is in progress. 



Acknowledgements We thank K. Balzer for useful 
discussions. 



Appendix A. Variational differentiation of 
abstract Hilbert space vectors 

In order to obtain the key equations for the UBHF 
approximation, one has to perform a variational deriva- 
tive with respect to all TV" orbitals. In this appendix 
we derive some differentiation rules that are needed to 
deduce (6). In the following the wavefunction of the n- 
th orbital will be denoted ifi n (x). Let us now consider 
the variational derivative of the matrix element of an 
arbitrary operator O: 
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S S f 

= <5in J 0{x,y)ipj(y)dy = Sin {x\6\j) 

= (-l(^(W>), (A.1) 

where the latter transformation can be regarded as the 
definition of the derivative with respect to an abstract 
Hilbert space vector: 



6(n\ 



(i\0\j) :=5 in O\j). 



(A.2) 



To verify that this definition makes sense we consider 
the same derivative but in an arbitrary representation. 
Therefore, we differentiate the given matrix element 
with respect to c* 7 - the 7-th expansion coefficient of 
the n-th orbital 



dc* 



(ipi\0\ipj) = 



d 

dc* 



a/3 



= E 



dc^ a 
3c* 



O a f3Cj/3 — S in ^2 7 pC 



■30 



<7|Ob) = <7l(^(i|Ob>)- (A.3) 



For the differentiation of the two-particle integrals, we 
can use the same ideas to derive the following equation 



S(n 



(ij\w\kl) \k}+5 jn J ik \l) (A.4) 

= S in K jk \l) + SjnKu \ k) . 



Thus we have some freedom to define the remaining 
operator. The first one that appears in this equation is 
defined above Eq. (8) whereas the exchange operator 
Kjk is given by 



Kjk = j 



dx dy tp*(x)w(x,y)(p k (y) \y) (x\ . (A.5) 



With these rules for the functional derivative of matrix 
elements, together with the product rule, 



5(n\ 



{Fi[- ■■ ,(n\ ,...]• Fz[. .. ,(n\ ,. . .]) 



\S(n\) 



F 2 + F-l 



SF 2 
S{n\ 



), (A.6) 



one derives equations (6) by differentiating the func- 
tional (3). 
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